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CERTIFIED REDUCED-BASIS SOLUTIONS OF VISCOUS BURGERS EQUATION 



PARAMETRIZED BY INITIAL AND BOUNDARY VALUES 
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Abstract. We present a reduced basis offline/online procedure for viscous Burgers initial bound- 
ary value problem, enabling efficient approximate computation of the solutions of this equation for 
' parametrized viscosity and initial and boundary value data. This procedure comes with a fast-evaluated 

■ rigorous error bound certifying the approximation procedure. Our numerical experiments show signif- 

icant computational savings, as well as efficiency of the error bound. 

1991 Mathematics Subject Classification. 35K20, 35K55, 65M15, 65M60. 



The dates will be set by the publisher. 



Introduction 



This paper is set in the context of sensitivity analysis and uncertainty analysis in geophysical models. Such 
models typically involve a wide range of parameters, such as: source terms (climatic forcings, heat/ wind/matter 
fluxes), boundary conditions (forcings, open boundaries), and the initial state of the system. 
| Their study generally leads to parametrized partial differential equations (PDEs) . These equations often involve 
t-H ■ poorly-known parameters. Therefore, it is important to be able to measure the impact of a given parameter on 
^ | the quality of the solution, and also to identify the "sensitive" parameters, that is, the parameters for which a 
• i-H ■ small variation implies a large variation of the model solution. Due to their ability to perform global sensitivity 
| analyses for nonlinear models, stochastic tools [13, 23] are rapidly expanding. These methods require "many 
queries," that is solving the parametrized PDE for a large (say, thousands) number of values of the parameters. 
When analytic solution to the PDE is not known (as it is often the case), one has to use a numerical method 
(such as finite difference or finite clement) to compute an approximate value of the solution. Such methods lead 
to computer codes that could take a large time to produce an accurate-enough approximation — for a single 
value of the parameter. Having the "many-query" problem in mind, it is crucial to design a procedure that 
solves the equation for several values of the parameter faster than the naive approach of calling the numerical 
code for each required instance of the parameter. 

The reduced basis (RB) method is such a procedure; we split the overall computation into two successive parts: 
one part, the offline phase, makes use of the standard, computationally intensive numerical procedure used to 
solve the PDE to gather "knowledge" about solutions of the latter; and the other one, the online phase, where 
we rely on data collected during the offline phase to compute, for each desired instance of the parameter, a 
good approximation of the solution, for a per-instance cost that is orders of magnitude smaller than the cost 
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of one run of the standard numerical code. The advantage is that, for a sufficiently large number of online 
evaluations, the fixed cost of the offline phase will be strongly dominated by the reduction in the marginal cost 
provided by the online procedure. This cost reduction is made possible by the fact that, in most cases, the 
desired solutions of the PDE, for all the considered values of the parameter, lie in some manifold of functions 
that is "close" to a low-dimensional linear subspace. One goal of the offline phase is to find such a suitable 
subspace, so that the online procedure can look for the solution of the PDE as an element of the subspace — 
so as to reduce the number of degrees of freedom and thus the computational cost. One interesting feature 
of the RB approach is that it comes with an online error bound, that is a (provably) certified, natural norm, 
fast-computed (i.e. almost of the same complexity of the online phase) upper bound of the distance between 
the solution provided by the online phase (called the reduced, or online) solution and the one given by the 
standard, expensive numerical procedure (called the full or reference solution). This "certified RB" framework 
has been developed for afjinely parametrized second-order elliptic linear PDEs in [19]. It has been extended to 
nonlinear, non-affinely parametrized, parabolic PDEs, see e.g. [11], [10] and applied to problems such as steady 
incompressible Navier-Stokes [25]. Moreover theoretical work has been done to ensure a priori convergence of 
the RB procedure [4]. 

In this paper, we are interested in the RB reduction of the time-dependent viscous Burgers equation (which 
will serve as a "test case" for the "real" equations modelling geophysical fluids we are interested in). Papers [12] 
and [22] extend certified RB methodology to linear initial-boundary value problems. The case of homogeneous 
Dirichlct boundary conditions, zero initial value and fixed (i.e., not parametrized) source term has been treated 
in [26], [18] ; in these works, the only parameter was the viscosity coefficient. Parametrization of initial and 
Dirichlct boundary conditions (treated using a conversion to a homogeneous Dirichlet problem) has been done in 
[16], for a general multidimensional quadratically nonlinear equation. Our methodology allows parametrization 
of the viscosity, the initial state, the source term and and of the boundary conditions. Compared to the works 
cited above, we use a weak (penalization) treatment of the Dirichlct boundary conditions. We will see that this 
weak treatment is more favorable in terms of both computation and storage complexities. Besides, our paper 
features a new error bound, which has shown to be, in all the testcases we performed, much more efficient than 
the existing bound. In particular, we will show that our bound exhibits improved sharpness for low viscosities, 
where [18] pointed out the moderate efficiency of the presented a posteriori bound. 

This paper is organised as follows: in the first part, we introduce the viscous Burgers equation, and present 
a standard numerical procedure used to solve it; in the second part, we expose our offline/online reduction 
procedure; in the third part, we develop a certified online error bound; finally in the fourth part we validate 
and discuss our results based on numerical experiments. 



In this section, we describe the model we are interested in. Subsection 1.1 introduces the viscous Burgers equa- 
tion, while Subsection 1.2 presents the "full" numerical procedure on which our reduction procedure, described 
in Section 2, relies on. 

1.1. Equation 

We are interested in u, function of space x £ [0; 1] and time t £ [0;T] (for T > 0), with regularity: 
u £ C 1 ([0, T], .ff^QO, 1[)), satisfying the viscous Burgers equation: 




where v £ R+ (R denotes the set of real numbers, the set of positive real numbers) is the viscosity, and 
/ £ C*° ([0,T],L 2 (]0,1[)) is the source term. 



1. Model 





(1.1) 
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For u to be well-defined, we also prescribe initial values uq G -ff^QO, 1[): 

u(t = 0,x)= u Q (x) Vx G [0; 1] 

and boundary values b ,bi G C ([0, T]): 

/«(*.» = o) = M*) WG[0;T] 

|n(t,x= 1) = 

Where feo 5 ^i and Mo are given functions, supposed to satisfy compatibility conditions: 

u o (0)=b o (0) and uo(l) = 6 X (0) 



(1.2) 



(1.3) 



(1.4) 



This problem can be analyzed by means of the Cole-Hopf substitution (see [14] for instance), which turns (1.1) 
into the heat equation, leading to an integral representation of u. 

1.2. Numerical resolution 

We now describe the "expensive" numerical resolution of the problem described above that will serve as our 
reference for the reduction procedure described in the next section. We proceed in two steps: space discretization 
in paragraph 1.2.1 and time discretization in paragraph 1.2.2. 

1.2.1. Space discretization 

For space discretization, we use a P 1 finite element procedure with weak (penalty) setting of the Dirichlet 
boundary conditions (1.3). 

We first have to write the weak formulation of our PDE ; to do so, we multiply (1.1) by a function v G -ff 1 (]0; 1[) 
and integrate over ]0; 1[: 



1 du 



(t, x)v{x)dx + 



1 f 1 d(u 2 ) 



dx 



(t, x)v(x)dx 



1 n2 



d u 
dx 2 



(t, x)v(x)dx 



f(t,x)v(x)dx W G H 1 (]0; 1[) Vi G [0;T] (1.5) 



Next, we integrate by parts the second and the third integral appearing in the left hand side of the previous 
equation: 



1 d(u 2 



d.v 

1 d 2 u 
o dx 2 



(t, x)v(x)dx 
(t, x)v(x)dx 



(t,x)^(x)dx+ [u 2 (t,-)v]l 



1 du dv 



! 1 



dx 



(t,-)w 



Inserting this into (1.5), we get: 

f 1 du . 1 f 1 2 dv f 1 du dv 

J 1u {t ' x)v{x)dX -2j + d- X {t > x) Tx {x)dX 



du 
dx 



(t,-)v 



= f{t,x)v(x)dx Vu G i? 1 (]0; 1 [) Vt G [0;T] (1.6) 
o Jo 



To get rid of the two boundary terms arising in the integrations by parts, one usually restricts v to satisfy 
v (0) = v(l) = so as to make the boundary terms disappear; the Dirichlet boundary conditions (1.3) are 
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then incorporated "outside" of the weak formulation. However, the reduction framework we are to expose later 
requires the boundary conditions to be ensured by the weak formulation itself. The Dirichlet penalty method, 
presented in [2], is a way of doing so, at the expense of a slight approximation error. This method entails 
replacement of boundary conditions (1.3) with the following conditions: 

1 du 

-~-u 2 (t, x = 0) + v—(t, x = 0) = P (u(t, x = 0)- b Q (t)) 

! 2 d t W£[0;T] (1.7) 

-u 2 (t, x=l)-v—{t,x = l)=P (u(t, x = 1) - 6i(*)) 

with a fixed penalization constant P > 0. 

The intuitive idea underlying (1.7) is that it can clearly be rewritten as: 



1 ( 1 du 
u(t, x = 0) = b o (t) + - i-^u 2 {t, x = 0) + i/^H(t, x = 0) 

u(t, x = l) = b 1 (t) + ^ (\u 2 (t, x = l)- i/^(t, x = 1) 



Vi G [0; T] 



so that (1.3) is asymptotically verified for P — > +oo. The reader can refer to [3] for rigorous a priori error 
estimates when using Dirichlet penalty in the linear elliptic case. 

In practice, we can check if our approximation is sufficiently accurate by means of the following a posteriori 
procedure: we take for P some large value (typically P = 10 ), we compute (numerically) an approximate 
solution Ud, using the procedure we are currently describing, and we check if an indicator of the amount of 
failure in verification of (1.3) is small enough; such an indicator can be, for instance: 

e b = sup[max(|u d (£,x = 0) - &o(*)|, \u d (t,x = 1) - h(t)\)] (1.8) 
t 

where the supremum is taken over all discrete time steps. If this indicator is larger than a prescribed tolerance, 
then P has to be increased. Our numerical results in Section 4 will assert this condition. We can then invoke 
the well-posedness of the boundary/initial value problem (specifically, continuous dependence on the boundary 
values) to ensure that the solution of (1.1), (1-2) and (1.7) will be close to the solution of (1.1), (1-2) and (1.3). 
This reasoning is analogous to the one made when omitting the approximation made when replacing exact 
boundary values by their discretized counterparts. 

Going back to our weak formulation, we multiply the first line of (1.7) by v (0), the second one by v(l) and add 
up these two equations. We get that: 



du , . 1 



P [(u(t, x = 0)v(0) - &o(*MO)) + (u(t, x = l)v(l) - &i(iMl))] 



Putting it back into (1.6) and isolating the terms not involving u on the right-hand side yields the following 
weak formulation: 

1 du 1 f 1 dv f 1 du dv , , „, , SN 

, Tt v -2] Q U dx + V L ^ + fWM = o)»(o)i,( ilI = iMD) 

1 

f(t,-)v + P(b o (t)v(0)+b 1 (t)v(l)) We-ff^lDVie [0;T] (1.9) 
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that is: 



( K — (t,-),v)+c(u{t,-),u(t,-),v) + va{u(t,-),v)+B(u(t,-),v) 

= t(v,t) + b (t)p (v)+b 1 [t)p 1 (v) V«Gff 1 aO;lD.Vt6[O j r| (1.10) 

by introducing the following notations (for all v,w,z £ -H^QO; 1[), t e [0; T}): 



and: 



B(w, v)=P (w(0)v(0) + w(l)v(l)) 
Po(v) = Pv(0) 

c(w, v, z) = 



i(v,t)= / f(t,-)v 
Jo 



dz 
dx 



7r : < 



for every w, v and z in -ff 1 (]0; 1[) for which wv— G i 1 (]0; 1[). 

ox 

The weak formulation is then discretized with Lagrange P 1 finite elements (see [21]) by choosing some integer 

Af and considering a uniform subdivision of [0; 1] with j\f + 1 nodes: {xi} i=0 ir and, for each i = 0, . . . ,Af we 

denote by </>,; the piecewise-affine "hat" function whose value is 1 on Xi and on every other nodes. 

We denote by X the linear subspace of f/^QO; 1[) spanned by {</>i} i=0 jy. We also set ||-|| to be the L 2 (]0; 1[) 

norm. 

Every ip € X can be written as tp = ipj&jt an d we have ip(xi) = ipi for all i = 0, . . . , TV. This justifies that 

3=0 

7r defined below is a projection of ff 1 (]0; 1[) onto X: 

H\]Q;l[)^X 

M 
3=0 

The space discretization of our problem is the following: for all t £ [0; T], find u(t, •) E X so that : 
u(t = 0, •) = tt(uo) 

(j£(t, •), «) + c(u(t, •), u(t, •), «) + w(u(t, •), «) + fl(u(t, ■). «) ( LU ) 

= 4(w,t) + 6o(*)A)(w) + iiWAM Vu e x 

where 

«*(«,*) = / T (/(<»•))« 

Jo 

Note that u now stands for a discrete solution, while it was used to designate an analytic solution before. 
1.2.2. Time discretization 

We now discretize (1.11) in time using the backward Euler scheme: we choose a timestep At > and consider 

T 

an uniform subdivision of [0;T]: {tk = fcAt} fc=0 j- where T = 
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Our fully discrete problem is: for k = 0, . . . , T, find u € X, approximation of u(tk, ■), satisfying: 

u° = n(u ) (1.12a) 

and: 

k k — l 

(~ — ~Ki ' W ) + C ^ fc ' ^ + va ( uk ' v ) + B (u k ,v) 

= lAv,tk) +b (t k )f3 (v) +b 1 (t k )p 1 (v) VveXVk = l,...,T (1.12b) 

We sequentially compute {u fe }fc=o,...,T m the following way: w° comes straightforwardly from (1.12a), and for 
fc = 1, . . . , T, w fc depends on it*^ 1 through (1.12b), which can be rewritten: 

^(u fc , w) + c(u k ,u k , v) + va(u k ,v) + S(u fe , «) 

= -^ t (u k - 1 ,v)+l w {v,t k ) + b (t k )p (v) + b 1 {t k )p 1 (v) Wei (1.13) 

We can now expand our unknown «' e I on the {<f>j}j basis: u k = u k fai, and the vector (u k ) . becomes our 

3=1 

new unknown. 

Moreover, it is sufficient for (lincar-in-u) relation (1.13) to be satisfied for all v in a basis of X, namely for 
v = fa, Vi = 0, . . . , J\f. 

So (1.13) can be rewritten as a nonlinear (due to the nonlinearity in c(u k , u k , v)) system of Af + 1 equations (one 
for each instantation v = 4>i) involving (v!f) ._ Q This nonlinear system is solved using Newton iterations: 

N 

starting with an initial guess u k , one looks for S = Sj(f>j so that u k = u k + S satisfies the linearization near 
8 = of (1.13) for v = 4>i, i = 0, . . . ,Af, that is to say: 

-£-(p}° + 6, 4>i) +c(t?,u*, ^<)+2c(t?, 6, fa) + ua{v* + 8, fa) + B(v}* + 5, fa) 

= ^(u k -\fa) + e w (fa, t k ) + b (t k )p (fa) + h(h)M<l>i) Vi = 0, . . . ,7V (1.14) 

System (1.14) is a (TV + 1) x (J\f + 1) linear system involving (Sj),_ Q Once solved for 8, one can test for 
convergence of the Newton iteration: if the norm of 8 is smaller than a prescribed precision, then we stop here, 
produce u k and get to the next time step k + 1; otherwise we do one more Newton step, this time using u k as 
initial guess u k . 

Wc can note that the linear system to be solved at each Newton step is not symmetric but is sparse. Due to this 
sparsity, its solution (using an iterative method such as BICGSTAB or GMRES) takes about O (M 2 ) operations 
in the worst case. One can take advantage of the tridiagonal structure of the matrix (which is present in the 
one-dimensional case, since fa and faj have no common support if \j — i\ > 1) and use Thomas' algorithm [24] 
to solve the system with 0(N) operations. 

2. Reduction procedure 

In this section, we show the offline/online procedure announced in the introduction to produce reduced basis 
solutions of the problem formed by (1.1), (1-2) and (1.3), based on the "full basis" numerical method presented 
above. We begin by a description of our parameters in Subsection 2.1. Our offline/online reduction procedure 
is described subsequently, in Subsection 2.2. 
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2.1. Parameters 

We parametrize Uq, bo, b\ and / as: 

n(bo) n(bi) 

b (t) = b 0m + E A bo <P bo (t) h(t) = b lm +Y, 4 I$ i 6l (*) 

;=i 2=1 

«t(/)"s(/) n(« ) 

= / m + e E AW T (*) $ f (*) = u °™+ E ^r o$ r°w 

i=X p=l (=1 

The functions involved in these linear combinations Q b ° , ^j? 1 , $f T , and tfr™ as well as the number of terms 
in the decompositions n(bo), n(bi), nr{f), ns{f) and n(iio) are fixed, while our parameters, namely: viscosity 

v, coefficients b 0m , b lm , f m and u 0m , and (A? ) , [A bl ) , (a£) and 

V 1 /J=l,...,n(6 ) V 1 /J=l,..., n (6i) V lp / 2=l,...,n T (/);p=l,...,n s (/) 

(A"°) J=1 , } live in some Cartesian product of intervals V' , subset of R 1 + 4 +»(M+MM+™t(/)« s (/)+«(«o)^ 
Note that m stands for mean and is not a "numerical" index. 

The compatibility condition (1.4) constraints bo m and b\ m as functions of the other parameters: 

n(«o) 

&Om=«0m aild ^lm = U 0m + ^ ^4"°$"°(x) 

{=1 

so that our "compliant" parameters actually belong to V defined by: 



1,"S(/)' 



— |m — { v i bomi Aj , ■ • ■ , A b ^ bg ^ bi m , A bl , . . . , A*° bi ^, f m , A^,A^2, ■ • ■ , A^ t 

4.1, • • • , 4,« s( /)> ■ • • - <(/),n s( /)^0 m , A?", ... , A^J € 7" satisfying (1.4)} (2.1) 



2.2. Offline/online procedure 

The key heuristic [19] for RB approximation of the (linear) parametrized variational problem is the following: 
given /i G V, 

find u(fx) € X so that A(u(fi), v; fi) — L(v; fx), Vu £ 1 
is to choose a parameter- independent family 1Z of linearly independent functions in X — with jj=TZ <C dimX, 
to achieve computational economy — and then, given an instance of the parameter, to search for the reduced 
solution: 

u{ijl) G SpanT?., so that A{u{fj), v\ /j) = L(v; /i) \/v G SpanTvL 
Let us apply this idea on problem (1.12). We rely on the procedure described in [18], modified to allow 
parametrization of initial condition, boundary data and source term. To simplify notations, we do not explicitly 
write dependence of u and u on \i. 

We suppose that a reduced basis TZ = • • • , Ov} has been chosen (see Subsection 2.3 for one way to do so); 
we define X = SpanT?. and we look for {u k } k _ T C X satisfying: 

u° = n(ir{u )) (2.2a) 

and: 

1 \ 

— ,vj + c(u k ,u k ,v) + va(u k ,v) +B(u k ,v) 

= £*(v,t h ) + b (t k )p (v)+b 1 (t k )p 1 (v) V«GXVfe = l,...,T (2.2b) 
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where tt is the orthogonal projection from X onto X, with respect to the standard L inner product on X 



(w, v) = wv. 
Jo 

The offline/onlh 

2.1: since the constant function 1 = <pj belongs to X, we have: 



The offline/online procedure for computation of u will come easily from our parametrization of uq in Subsection 

TV 



3=0 

n(uo) 



U° = U Qt - 



,7r(l) + E ^"^WO)) (2-3) 
(=1 

We now discuss computation of u k from u fe_1 for k = 1, . . . ,7". We are willing to proceed with Newton steps 
as for the solution of (1.12b). We denote, for k = 1, . . . ,T, by 

3=1 3 = 1 

respectively, the reduced solution at time and previous guess for the reduced solution at time t^. 

Our procedure relies on the following proposition: 

N 

Proposition 2.1. (1) The Newton increment 8 = V^<5jCj satisfies the following equations: 

3=1 

E *i I + 2 E 4 c (0'. c-, co + mq,q) + co 



At 

3=1 I 3'=1 



= E w 3 _lA! 7^ + MCi, *fe) + &o(*)A>(C0 + 6i(*)^i(Ci) 



3=1 



(2) We /icrae: 



S u 3 fc (%r- + E4 c (C3'>C3,co+^(C3,co+s(Ci,co ] v* = i,...,jv (2.4) 

3=1 \ 3"' = 1 



»1 «t(/)"s(/) .1 

4(Ci,* fe ) = /m / 0+ E E (**) / *(K S (•)) c< (2.5) 
- /o ;=i p=i J ° 



/or all i = 1, . . . , N and k = 1, . . . , T. 
Proof. (1) Equation (2.2b) for w fe = M fc + <5 linearized near 5 = 0, for v = £i, Vi = 1, . . . ,N is: 

+ (5, Ci) + c(5*, CO + 2c (^> <*> CO + ^«(^ + CO + + 5, CO 

= ^<«* -1 jC0 + 4(Ci.*fc) + 6o(*)A)(C0 + &i(*)/MC0 Vi = 1, • • • ,N 

Rewriting this equation using expansions of u k and 8 in 1Z and linearity of (•,■), c, a and B with respect 
to their first argument, and putting all (5j)j -dependent terms on the left-hand side give the announced 
equation. 
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(2) is a direct consequence of the parametrization of / given in Subsection 2.1. 

□ 

The following Proposition 2.2 justifies the wcll-dcfincdness of the Newton iteration, for an orthonormal reduced 
basis . . . , £n}- In practice, the Gram-Schmidt process can always be used to ensure this condition. 

Proposition 2.2. Suppose that {£i, . . . , £jv} is orthonormal with respect to (•, •). 

Then, for At small enough, i.e. At < At* (i>, {Ci}i=i,...,iv)j and initial guess u k sufficiently close to u k , the 
Newton iteration (2.4) is well defined and converges (quadratically) to u k . 

Proof. Iteration (2.4) is a Newton iteration for solving F(x) = a, for appropriate a and F given by: 

F(x) = J2 + c (*> ^ O + M*, Q) + B(x, 0)) & 

i— 1 ^ ^ 

We apply the result stated pp. 353-355 of [21] and pp. 362-367 of [6]. To do so, we have to check that, for At 
small enough: 

(1) the differential of F at u k , denoted by DF(u k ), is invertiblc; 

(2) DF(-) is Lipschitz-continuous i.e. there exist L > so that 

Vi.i'el, Vuel, \\DF(x) ■ v - DF{x') ■ v\\ <L\\x-x'\\\\v\\ 

It is easy to check that the matrix of DF(u k ) in the reduced basis . . . , (at} is diagonally dominant, hence 
invertible, for: 

At < At := mm — r , 

4=1 '-' JV £f =1 I2c(w fc , 0, CO + va(Q,(i) + B(Q,Ci)\ 

and (2) is a straightforward computation. □ 

We put our offline/online procedure in Algorithm 1. Let us make some remarks about the complexity of the 
above online algorithm, in contrast with the "full basis" one described in Section 1.2: 

Remark 2.3. (1) The most computationally demanding step is (2) (c), since it involves resolution of a 
(nonsymmetric, dense) N x N linear system; the "full basis" counterpart solves (AT + 1) x (J\f + 1) 
nonsymmctric tridiagonal system. Thus significant computational savings are expected for N <C M . 

(2) Thanks to our parametrization of /(t, •), all integrals over [0; 1] in equation (2.4) can be precomputcd 
during the offline phase, yielding a Af- independent online phase. This means that one can in principle 
choose arbitrary high precision on the full model without impact on the marginal cost of evaluation 
of an online solution. This "A/'-independence" property shall be required of every complexity of any 
online procedure. One should also note that our online procedure does not produce "nodal" values 
u k (xi), i = 0, . . . , Af (as this would clearly violate the A/'-independence), but rather the components of 
u k in the reduced basis. 

(3) Taken independently, the number of parameters n(uo), n(bo), n(bi), n$(f) and nx(f), as well as the 
number of timesteps T have a linear impact on the online complexity. Moreover, due to the double sum 
in (2.5), the online complexity is proportional to «s(/)«t(/). We note that an advantage of treating 
the Dirichlet boundary condition weakly, as we do in this paper, is that nr{f) does not get increased 
by functions of n(bo) or n(bi). This is a clear advantage of our method: when boundary conditions are 
treated by returning to an homogeneous Dirichlet problem, as in [16], n(bo) + n(b\) terms are added in 
the parametrization of /. 
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• offline: 

(1) choose a parameter-, and time-independent reduced basis {£1, . . . , (n} (see Section 2.3) 

(2) compute and store the following parameter-independent functions of X and scalars, for all i, j, f 
l,...,N,l = l,..., n(u ), p=l,..., n s (f): 

tF(I) 7r(7r($r°(0)) 

<G:0> a(C i; Q) 

c (C,v,0, 0) 

A,(0) /3i(0) 

ni o ok, 

o Jo 

• online: 

(1) assemble u° as the linear combination (2.3) ; 

(2) forfc = l,...,T: 

N 

(a) set up an initial guess u k = u jCj > 

i=i 

(b) compute and store 4r(0>*fc) by using (2.5) ; 

JV 

(c) look for S = 8jQ by solving the linear system (2.4) ; 

_ 3=1 

(d) set w fc w fc + J ; 

(e) if || 5 1 is small enough: 

(i) output u 

(ii) set k = k + 1 

(f) else: 

(i) update the guess : u k <— u k , i.e. u k j «— Vj = 1, . . . , N 

(ii) go back to (c) 

Algorithm 1: offline/online procedure 



2.3. Choice of the reduced basis 

In this subsection, we describe different ways of choosing a pertinent reduced basis {Ci, ■ ■ • , 0v}- These lead to 
three different bases fitting into the certified (that is to say, the three admit the same procedure for online error 
bound we describe in Section 3) reduced basis framework. 

The first is based on proper orthogonal decomposition (POD), the second is based on a "greedy" selection 
algorithm. The third is an hybridation of POD and greedy. These methods are standard in the literature. 



Notation 

The two procedures described below involve computation of the reference solution for different instances of the 
parameter, so we should use special notations, local to this section, to emphasize the dependence of the reference 
solution on the parameters. We define a parametrized solution u by: 

({l,...,T}xP^X 
u : < , 

(k,n) i — } u(k,/i) — u satisfying (1.12b) for \x as parameter 
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2.3.1. POD-driven procedure 

We denote by {uj(n)}j the coordinates of u(k, /u) in the basis {<f>o, . . . , 4>Af} : 

N 

u{k,n) = ^2 l u k j {ix)4> ] 

3 = 1 

In the POD-based procedure (see [5]) of the reduced basis choice, we choose a finite-sized parameter sample 
5 C V, compute the reference solutions u(k,/j,) for all k = 1, . . . ,T and all fi € 5, and form the snapshots 
matrix containing the components of these solutions in our basis {4>o, ■ ■ ■ > 4>J\f} '■ 



M 



( «o(Ml) w o(Mi) 



K (Ml) «o(M2) 



wo (ms)\ 



u if(ps)J 



where S = ■ ■ ■ ,Ms}- 

One can check that M has A/" + 1 rows and S(T + 1) columns. 

To finish, we choose the size N < S(T + 1) of the desired reduced basis, we form the S(T + 1) X S(T + 1) 
non-negative symmetric matrix M CIM, where f2 is the matrix of our inner product (•, •), to find z\, . . . , zjv 
the N leading nonzero eigenvectors of this matrix (that is, the ones associated with the N largest eigenvalues, 
counting repeatedly possible nonsimple eigenvalues), and, for i = 1, . . . , N, the coordinates of Q with respect 
to {</) , ... , 4>at} are given by: 

uh\ Ml ' (2 - 6) 

2.3.2. "Local" greedy selection procedure 

The greedy basis selection algorithm (cf. [9,11]) is the following: 



Parameter: N, the desired size of the reduced basis. 

(1) Choose a finite-sized, random, large sample of parameters ScP. 

(2) Choose \i\ £ V and k\ & {0, . . . , T} at random, and set 



Ci 



u(Hi, fci) 

|ti(/il,fcl)|| 



(3) Repeat, for n from 2 to N: 
(a) Find: 



(Pn,k n )= argmax e*{p,k) 

(A«,fc)eHx{0,...,T} 



where e*(/i, k) is a (fastly evaluated) estimator of the RB error ||m(/j., k) — u(fi, k)\\, where u(/i, k) 
stands for the RB approximation to it(/x, k) using {Ci, . . . , Cn-i} as reduced basis (see below). 

(b) Compute (* = u(p n , k n ). 

(c) Using one step of the (stabilized) Gram-Schmidt algorithm, find €E Span{^x, • • • , Cn-i> Cn} so 
that {Ci, . . . , Cn-i, C«} is an orthonormal family of (X, (■, ■)). 



Algorithm 2: Greedy basis selection 

The "greedy" name for this algorithm comes from the fact that the algorithm chooses, at each step of the repeat 
loop, the "best possible" time and parameter tuple to the reduced basis, that is the one for which the RB 
approximation error is estimated to be the worst. 
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Let's now discuss the choice for the error indicator e* . A natural candidate would be the online error bound e 
described in Section 3. However, as we shall see in the next section, the bound for timcstep tk is a compound 
of the "propagation" of the error made in the previous timesteps tojti, . . . ,tk—i (the Ek-i term) and the "local 
error" just introduced at the fc-th time discretization. Hence, this error estimator has a natural tendency 
to grow (exponentially) with k. Thus using it as error indicator e* will favor times k n near final time T to 
be chosen at step 3. (a) of the algorithm below. Such choices are suboptimal, because including them in the 
reduced basis will not fix this exponential growth problem which is inherent to our approximation. Instead we 
use, as in [11], a purely local-at-tk error indicator, that is: the computable error bound described in Section 3 
when taking £k-i = 0. 

It has been noted in the literature that the greedy procedure can "stall" (i.e. select vectors for addition in the 
basis which have no effect in decreasing the error bounds). The author in [9] provides a "back-up procedure" in 
such a case. In our experiments, however, we have not noticed any "stalling" of the greedy procedure, maybe 
because of our sharper error bound. 

2.3.3. POD-Greedy procedure 

We can also make use of the POD-Greedy procedure [12] in Algorithm 3. This procedure aims to combine 
the advantages of the greedy and the POD based procedure; it has also been proposed so as to overcome the 
"stalling" of the greedy procedure. Again, an error indicator e* has to be chosen. We do as in [17], where the 



Parameter: P\. 

(1) Choose a finite-sized, random, large sample of parameters ScP. 

(2) Choose Hi £ V at random, and choose as the current reduced basis B an orthonormalized basis of 

Span{it(/ii,io),ti(Mi,*i), ■ ■ -,u(hi,T)}. 

(3) Repeat, until the desired number of items in the basis is reached: 

(a) Find: 

/a* = argmaxe*(/i) 

where e*(/i) is a (fastly evaluated) estimator of the RB error \\u(fi) — u((i)\\, where u(/j,) stands for 
the RB approximation to u(jx) using B as reduced basis. 

(b) Append to B the Pi leading POD modes of {u proj (/i*, to), w proj (/A fi), . . . , u proj (^*, T)}, where the 
proj superscript denotes projection on the orthogonal complement of Span(S). 

(4) Output B = {d, . . . , C##} as reduced basis. 

Algorithm 3: POD-Greedy basis selection 



authors take the online error bound at final time. 
2.3.4. Expansion of the basis by initial data modes 

In case they did not get chosen by the POD or greedy algorithm, a classical strategy is to initialize the basis 
selection algorithms by taking in the reduced basis the constant function 1, and the functions 3>"°(-) for I = 
1, . . . , n(ua). This may increase the size of the reduced basis (and thus online computation times) but ensures 
that u° = uq (i.e. initial error is zero). Such an enrichment can possibly be a good move, as the error gets 
accumulated and amplified throughout the time iterations, zero initial error will certainly reduce the (estimated, 
as well as actual) RB approximation error. 

3. Error bound 

In this section, we derive a parameter and time dependent online error bound e h (for k = 1, . . . , T) satisfying: 
\\ u e ~ uk \\ — £k i w h crc I = ( / v 2 ) , and u\ is our reference "truth" solution in X satisfying the fully 
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discretized PDE with strong Dirichlet enforcement, that is: 



,k-l 



At 



«S(0) = b a (t k ) 
u k e (l) = h{t k ) 



(3.1) 



where X is the "homogeneous" subspace of X: 



X = {v e X st. w(0) = w(l) = 0}. 



Our error bound should be precise enough (i.e., not overestimating the actual error «* — u fc || too much, and 
approaching zero as N increases) and online-efficient (that is, admit an offline/online computation procedure 
with an A/"-indcpendcnt online complexity). 

We notice that our error bound measures the error between the reduced and the reference solution; it does 
not reflect the discretization error made when replacing the actual analytical solution of the Burgers equation 
with its numerical approximation (3.1). This is consistent with the fact that RB methods relies strongly on the 
existence of a high-fidelity numerical approximation of the analytical solution by a discrete one, hence regarded 
as "truth". 

Subsection 3.1 deals with the derivation of the error bound; this error involves quantities whose computation is 
detailed in Subsections 3.2 and 3.3. 



3.1. Error bound 

The sketch of the derivation of our error bound is the following: we first give a "theoretical" error bound 
in Theorem 3.1; we then replace the uncomputablc quantities appearing in this bound by their computable 
surrogates in paragraph 3.1.2. 

3.1.1. Theoretical error bound 

Notation. We suppose that the convergence tests appearing in the Newton iterations performed in Section 
1.2.2 and Section 2.2 are sufficiently demanding so as to neglect the errors due to the iterative solution of the 
nonlinear systems (1.13) and (2.2b). 

We now set up some notations : first the error at time tk- 



u\ -u k if k > 

7r(uo) — tt (n(u )) if k = 



the residual form r k , for v £ X : 



r k (v) =tn(v,t k ) - ^> fe -u k \v) - c(u k ,u k ,v) - va{u k ,v) (3.2) 



and the "Xo-norm" of the residual: 



kfcllo = SU P r k{v) 

vex ,\\v\\=i 



We introduce : 

ipkip, w) — 2c(u k , v, id) + va(v, w) 

and the so-called stability constants: 



C k = inf ipk(v,v) (3.3) 
t;ex ,|H|=i 
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To finish, we define: 

% = MO)|||0o|| + Ml)|||<MI ; Ok=2\C k \i) h ; £= sup v(±- 

vex , \\v ||=i V-/V 

(£ is finite because Xq is finite dimensional), and, finally: 

f k = S ( |e fc (0)| |^(0o, 0i)| + |e fc (l)| |Vfe(0AT, ) 



£ 2 (|e fc (0) | + |e fc (l)|) B 5 , 2 2 v M0)| 3 + K(1)| 3 
3 ; ^=3^(efc(0) +e fe (l) ) ; = 3 



0? = 3 ; Cife = q^(efe(0) + e fe (l) ) ; £/ 

We will also make use of the standard notation 



[x]_ = max(— x, 0) Va; € R 

The theoretical foundation for our error bound is the following theorem. As the technique developed in [1 8] did 
not fit our problem (because of the weak Dirichlet treatment, whose advantage has been shown in Remark 2.3 
point (3)), we developed a new strategy for obtaining this error bound. 

Theorem 3.1. If: 

^ + C fc -#>0Vfc = l,...,T (3.4) 
then the norm of the error \\e k \\ satisfies: 



( Bk 2A * 



ifV k >Q 
ifV k <0 



with: 



A 1 , r M . „ 2r; fc + |le fc _ 1 ||+g(0o,0 1 )(|e fc (O)| + | efc (l)|) , B 
At ^ fe ' At ^°k + fk + \\r k \\o + £k 



rjk \\e k -i\\ +g?7fc((/)o,0 1 )(|efc(O)| + |e fc (l)|) 2 , 1 1 m 3 ^31 , C 7 , 11 11 

7fc = ^ + Vkfk + [C k \_r, k + - |e fc (l) -e fc (0) | + £ fe ' + ||r fc || % 

and: 

V k = (B k f+4A k ~f k . 
The proof of Theorem 3.1 is presented in the appendix. 
3.1.2. Computable error bound. 

We now find an efficiently computable (that is, with an offline/online decomposition, with a complexity of the 
online part independent of Af) error bound e k derived from the one described above; to do so we discuss each 
of the ingredients appearing in its expression. 

• Computation of the norm of the initial error ||eo|j is addressed in the next Section 3.2; the one of ||t>||q 
is in Section 3.3. 

• We have, for w G {0, 1}: 

e k {w) = u k -(w) ~ u k (w) = b w {t k ) - u k (w) 
so that efe(0) and efc(l) can be computed during the online phase. 
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• Similarly, the scalars ||0o||> ipki&Oi&i) an d ipk{4>N,N-i) can straightforwardly be computed online. 

• The "continuity constant" £ can be computed offline and stored by solving the optimization problem 
defining it. Thus rjk, and £ k can be computed online. 

• The exact value of Ck could be found by solving a generalized eigenvalue problem on X: Ck is the 
smallest A € R so that there exists z £ X , \\z\\ = 1 satisfying: 



i>l ym {z,v) = \{z,v) Vv£X 



with ym the symmetric bilinear form defined by: 



ip k vm (w, v) = va(w, v) + c(it fc , w, v) + c(ii fc ,v, w) Ww, v £ Xq 



The cost of doing so is prohibitive as it is an increasing function of dimX = J\f + 1. Instead, we will 



see in Section 3.4 how to compute lower and upper bounds C l k ^ and C k 
We can then compute the following lower and upper bounds for Ak- 



sup , j^tinf 



or < c k < c k 



sup 



A7 f = ^ + c k 



1 



inf 



1 



a sup - | ^sup _ £ 

"** ~ At + k 4 



and the hypothesis (3.4) is ensured by checking that A^ > 0. 
• We can also compute an upper bound of a k '- 



sup 



2?7 fe max(|C i 



sup I 

k 



• To compute an upper bound of Bk, we need to replace the preceding error norm || e*; i |j which is (except 

for k = 1) not exactly computable, with the online upper bound £k-i > ||efe_i|| at the preceding time 
step: 

2% + e fe _i 



B sup 



g(0o, 0i)(|e fc (O)| + |e fc (l)|) sup , .. .. B 

r a k + Jk + ll r fc|lo + Sfc 



At 



And 7fc gets replaced by its upper bound 7 fc 



sup . 



sup VkEk-i + <?u)(|efc(0)| + |e fc (l)|) 
7 fc * = ^ + »?*/«, + 



(jinf 



V 2 k + Q\e k (l) 3 -e k (0f\+C k + \\r k \\ Q Vk 



• We finally compute an upper bound for T>^\ 



q-ySUp 



(£f p f + 4A k lf j s k up if 7 f p < 0. 



and our "computable" error bound is then: 



( Jgsup 



S U /.) 

k 



B s k up 
A in f 



2A k nf 



if V s k up > 



if V s k up < 0. 



The remainder of the section consists in the description of computation of the four left-out quantities ||eo||, 
\\r k \\ Q ,Cr andC* 1 '. 
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3.2. Initial error 

The present subsection deals with efficient computation of the ||eo|| term in the computable error bound described 
in paragraph 3.1.2. We denote by H the Gram matrix of the family: 

{l - 5r(l), tt (*?(.)) - 5F (* ($?(.))) (C (tto) (-)) - 5? (vr (*^,(-))) } 

that is, H is the (1 + n(uo)) x (1 + n(uo)) symmetric matrix of all the inner products between two any of the 
above vectors and by eo the vector containing the components of eo with respect to the family above, i.e.: 

e =(u 0m ,Ar,...,A% Uo) ) T 

We have: 

Lemma 3.2. The norm of the initial error ||eo|| is given by: 

||eo|| = v / eo T He (3.5) 

Proof. Paramctrization 2.1 gives: 

n(uo) 

e = n(u ) n (n(u )) = u Qm (1 - ?r(l)) + £ A T (* W (")) " * W (')))) 

and the result follows from the expansion of ||eo|| 2 when eo is replaced by the expression above. □ 

This formula allows us to compute the time and parameter-independent Gram matrix H during the offline 
phase, and, during the online phase, to assemble the (1 + n(M ))-vector e and to perform (3.5) to get ||eo|| 
with an online cost dependent only of n(uo). 

3.3. Norm of the residual 

We now present the computation of the ||rfc|| term in the computable error bound of Section 3.1.2. 

• Let the expansions of the reduced solutions with respect to the reduced basis be: 

N 

u p = J2 uP jQ forpe {fc- l,fc} 

• Let G be the (1 + n s (f) + 2N + N 2 ) x (1 + n s {f) + 2N + iV 2 )-sizcd Gram matrix of 

/-pint p/S p/S p() -p{) -pc -pc pc pc pc pc -pa pa \ 

I 1 ' 1 j • • • ) 1 ns (/)i 1 1 ! • • • ) 1 JV' 1 1,1' 1 l,2i • • • J - 1 - l,iV) - 1 - 2,1: • • • ) 1 2,7V > ■ ■ • i 1 TV, AM 1 1 1 ■ ■ ■ ; 1 N S 
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where T int , T^ s , vf , T^j, , T™ e X (p = 1, . . . , n 8 (f)]j,j' = 1, . • • , N) satisfy: 

(r mt ,v) = f v Vvex 
i 



o 

(rf,v) = (Q,v) VveX 

(Tl f ,v) = c(Q,Q>,v) Vt-eXo 

(T^v}=a(Q,v) VveX 

(Those T's exist by virtue of the Riesz representation theorem). 
• Let p k be the following vector: 



™t(/) 

Pk 



f m , A{ j ${ T (t k )(j = l,...,ns(f)),-^(u>]-u>;- 1 )(j = l,...,N), 



i=i 



Lemma 3.3. We have: 



u^uf (j,l = l,...,N),^ (j = l,...,N))' 



IWIo - \/PfcGp fc (3.6) 
Proof. From the Riesz representation theorem, there exists a unique p k € Xq so that 

(p k ,v) =r k (v) VveX (3.7) 

and we have: ||r fc || = \\p k \\. 

From (3.7) and the definition of r k (3.2), we have that p k is defined uniquely by: 

(pk,v) = lx[v,t k ) - -^(u k - u k ~ x ,v) -c(u k ,u k ,v) - va(u k ,v)Vv e X 
because fto(v) = f3i(v) = B(-, v) = for all v s Xq. 

Using parametrization (2.5) of £„■(•, t&), we get that (3.7) is equivalent to: 



M = f m v+j2 E<*r(**)/ ^(^ s (-))^- X7 E(4-^ 1 )(o^) 

1=1 p=l J « A S = i 



AT /TV 

E u i E 4 c «j > G' > w ) + ^ . «) ] W e (3.8) 

3=1 \3'=1 



By the superposition principle, p k can be written as the linear combination: 



N N N 



nsU) /«!•(/) \ 1 JV 

. fc = / m r- + E E 4.*f (**) H 8 - ^E("3 - »Dr? - E E <4i ;,' -E w i r " 
p=i \ (=i y j=i 



3=1 3'=1 3=1 
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Thus p k contains the components of p k with respect to the family whose G is the Gram matrix, and so: 



\r k \\ =\\ Pk \\ = JplGp 



(3.9) 

□ 



The offline/online decomposition for computation of ||7"fc|| is as follows: in the offline phase, we compute the T's 
vectors, and compute and store their Gram matrix G. In the online phase, we compute p k and compute ||?~fc|lo 
using (3.6). Note that one can reduce offline and online computational burden, as well as storage requirement, 
by noticing that r| = 1^ • for all j,f. 

The cost of computation (and storage) of p k and \\pk\\ asymptotically dominates the cost of the online phase. 
This cost is in O ({ns(f)nr(f) + iV 2 ) 2 ^. Again, we see that our weak Dirichlet treatment, as nr(f) remains 
independent of n(bo) and n(bi), allows for a better complexity of the online phase. 

3.4. Lower and upper bounds on stability constant 

To find lower and upper bounds on C k efficiently in order to use them in our computable error bound of Section 
3.1.2, we turn to the successive constraints method (SCM) [15,18]. Here we present the application to our case 
for the sake of self-containedness. Our difference is the use of the metric (3.10) during the constraint-selection 
phase. 

Notation. As in Section 2.3, we will need to handle reduced solutions of several values of the parameter tuple 
p S V (see (2.1)), for different tiniest eps k = 1, . . . ,T. We thus define an application that is the "reduced" 
counterpart of u defined in Section 2.3: 

_ ( {l,...,T}xV^X 

I (k, p) i y u(k, p) = u satisfying (2.2b) for p as parameter 

We make C k depend explicitly on p by defining: 



C k (p) 



inf [2c(u(k, fi),v, v) + v(p)a{v : v)} 
vex , \\v =i 



SCM lower bound. We now proceed to the derivation of the SCM lower bound of C\(p). We use the RB 



N 



expansion: u{k,p) = Uj(p)Q to rewrite C k (p) as: 



3 = 1 



C k (p) = inf 

eex ,||«||=i 



JV 



^22uj(p)c{Q,v,v) + va{v,v) 



inf 

y=(yi,---,VN+i)ey 
inf J (p,k,y) 



N 



where: 

y = {y = (yx,...,y N+ i) e n N+1 \3v e x ,\\v\\ = l s.t. y 3 = c(Cj,v,v)\/j = 1,... 

and: 

N 

J{p,k, y) = 2 ^ u) (p) yj + vy N+1 

3=1 

We define, for a given "constraint subset" C C {1 T} x P that will he chosen later: 



,N,y N+1 = a(v,v)} 
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y=^(yi,...,y N+ i)e J[ [of* 1 ;^] W, k',y) > C k , (,/), V(//, A;') e <S(//, fc)| with<S(/i, A;) stand- 
ing for the set of the M points in C that are closest to (/x, k) with respect to this metric: 

dimP / i li \ 2 /U U>\ 2 

d{( Ml kW,k>))=Y, / _^ + V 1 (3 - 10) 

i= l \f J 'min f^max ) \ 1 J 

where (fx 1 , ... , /x dlmP ) are the coordinates of [i 6 V, and: 



Mmin = min /x , /i maK = max /! 

for i — 1, . . . , dim? 3 (here dimV = 1 + 2 + n(6o) + n(bi) + nT(f)ns(f) + n(uo) is the number of 
parameters). 

The metric defined in (3.10) quantifies proximity of two parameter-time tuples, with appropriate 
weighting so as to account for scaling differences between the parameters. 
• We further define: 

of" = inf c«i, v,v), Vi = 1, . . . , N o^S = inf a(«, v) 

vex ,\\v\\=i ueX ,IMI=i 

0™" = sup c(Ci,v,v), Vi = l,...,JV o^= sup o(«,t>) 

ueX ,IMI=i «ex ,||u||=i 

The SCM lower bound is then given by the following lemma: 

Lemma 3.4 (Proposition 1 in [15]). For every C C {1, . . . , T} X P and M € N, and every k = 1, . . . ,T, a 
lower bound for C k {n) is given by: 

Ci nf (ri=M[J(»,k,y)} (3.11) 

An algorithm for choosing a constraint subset C will be given after the description of the SCM upper bound 
and the SCM offline/online procedure. 
SCM upper bound. We define: 

y up = {y*(ki,fii) : i 1 /: (ki,fii)€C} 

where: 

C = {(fci,/ii), (fc 2 ,/i 2 ), ■ ■ • , 

and: 

y*{h,pLi) = arginf [J{(m, h,y)} (i = l,...,I) 
yey 

Lemma 3.5 (Proposition 1 in [15]). For every k = 1, . . . , T, an upper bound for Ck is given by: 

Cr»= inf J(^Ky) (3.12) 

yeyup 

SCM offline/online procedure. Relying on Lemmas 3.4 and 3.5, our offline/online procedure for computing 
C\ n and C s k uv reads: In the lower bound online phase, the optimization problem required to solve is a linear 
programming problem (LP) with N + 1 variables and N + 1 + M constraints (M one-sided inequalities and 
N + 1 two-sided). There are algorithms, such as the simplex algorithm (see [20] for instance), which solve such 
optimization problems under (on average) polynomial complexity with respect to the number of variables and 
number of constraints, even if they can be exponential in the worst cases. What matters here is this complexity 
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• offline: 

(1) choose M and constraint set C (see next paragraph) ; 

(2) compute and store cr™" and a™ ax (i = 1, . . . , JV+ 1) by solving a generalized cigenproblem on Xq 

(3) for each {k',fx') G C: 

(a) solve a generalized eigenproblem on Xo to find Ck'(fjf) (and store it); 

(b) let w G Xq be a unit eigenvector of the above eigenproblem; compute and store (in y up ) 
y*(k',fj,') using: 

y*{k',n')j =c((j,w,w) (j = l,...,N) 
y*(k' ,^') N +\ = a(w,w) 

• online: 

— for the lower bound: 

(1) assemble and solve optimization problem (3.11) ; 

— for the upper bound: test one-by-one each element of y up to solve (3.12). 

Algorithm 4: SCM offline/online 

is independent of J\f. The upper bound online phase has a complexity depending lincarily on the cardinality of 
the reasonably-sized C and on N. 

"Greedy" constraint set selection. To choose C in Algorithm 4, step 1, we can use the greedy constraint 
set selection Algorithm 5. The repeat loop in this algorithm can be stopped either when #C has reached a 



(1) choose M G N ; 




(2) initialize C = {(fci, Mi)} with arbitrary k\ G {1, . . . ,T} and \i\ G V\ 




(3) choose a rather large, finite-sized sample S C {1, . . . , T} X V ; 




(4) repeat: 




• using the "current" C to compute C sup and C m * , append: 




exp(C^( / ,))-exp(cr / ( / i)) 




(k*,n*) = argmax 

(*, M )es exp(C* fc "(/i)) 


to C. 



Algorithm 5: Greedy constraint set selection 



maximal value, or when the "relative exponential sharpness" indicator: 

exp(Cr^))-cxp(cfV)) 
(fe,,u)eH exp (C s k up {n)) 

gets less than a desired precision. We use a measure of the difference between the exponentials so as to account 
for the "exponential" effect of the stability constants on the error bounds [18]. 

As in the greedy algorithm for basis selection described at Section 2.3.2, this algorithm makes, at each step, the 
"best possible" choice, that is the value of the parameter and time for which the bounds computed using the 
current constraint set are the less sharp. 

A last remark we can give on the algorithm is about the trade-off in the choice of M: whatever M is, we always 
get a certified bound on Cfe(/i), but increasing M will improve sharpness of this bound, at the expense of an 
increase in online computation time. 
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4. Numerical results 

We now present some numerical results obtained with the methodology described above. We implemented it in 
a software package written in C++, using GNU OpenMP [8] as threading library, ARPACK [1] for eigenvalues 
computation and GLPK [7] as linear programming problems solver. 

For all the experiments above, the convergence test for Newton iterations when solving (1.13) and (2.2b) was 
the following: ||5|| < 3 x 10~ 16 . The penalization constant used was P = 10'. 
We also took M = #C = 10 as parameters for the SCM procedure. 

The $ functions appearing in the parametrizations of uo, bo, b\ and / are chosen to be sine functions with fixed 
known angular velocity. More specifically, we suppose that: 

n(6 ) n(i>i) 

b (t) = b 0m + A \ 8in(wf°t) h(t) = bim + A f M^t) 

i=i i=i 

»t(/)is(/) n(u ) 

f(t, x) = f m + ]T A l s H"{ T t) sin(c^' s x) uo(ar) = u 0m + £ A? sin(^ n x) 

i=i P =i i=i 

4.1. Reference solutions 

Figure 1 shows an example of the reference solution of (1.1), (1.2) and (1.7) every 10 timcstcps. The parameters 
were: 

Af = 40 At = .02 
T = 2 v = \ 

(4 1) 

b (t) = 1 h[t) = 1.28224 v ■ ' 

f(t, x) = l u„{x) = 1 + 2 sin(3a;) 
Figure 2 does the same with a lower viscosity. The parameters were: 



(4.2) 



N = 


40 


At = 


.002 


T = 


2 


v = 


.1 


bo(t) - 


1 


h(t) = 


1.28224 


f(t t x) = 


1 


u (x) = 


1 + 2 sin(3 



Figures 1 and 2 show the solution u of the viscous Burgers' equation plotted as functions of space x, for various 
times t, respectively for the parameter set (4.1) and (4.2). 

We have checked that the a posteriori indicator of boundary error (1.8) gets no higher than 6 x 10~ 7 in both 
cases. 

4.2. Reduced solutions 

Computational economy 

To show the substantial time savings provided by the reduced basis approximation, we compute the reduced 
solution for the parameters set given by (4.1), with Af = 60. The full solution (with Af = 60) takes 0.26s CPU 
time to be produced (when using Thomas' algorithm for tridiagonal matrices inversion). 

We use the POD-driven basis selection procedure to select the N = 7 leading POD modes (using S = 30 
snapshots). The resulting basis (with the functions sorted by increasing magnitude of eigenvalues) is shown in 
Figure 3. We did not make use of the "expansion" procedure described in Section 2.3.4. The overall CPU time 
for the offline phase was 6.36s. 

We used fixed parameters n(bo) = n(&i) = n${f) = nr{f) = n{v,o) — 1 and parameter ranges as shown in Table 
1. 
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0.2 0.4 0.6 0.8 1 

Figure 1. Full solution with high viscosity v = 1. Plots of the solution u of equation (1.1), 
(1.2), (1.7) as a function of space x, for various times t ranging from t = to t = 2 (the bottom 
lines correspond to high times). We use the parameters defined in (4.1). 




0.2 0.4 0.6 0.8 1 

Figure 2. Full solution with low viscosity v = .1. Plots of the solution u of (1.1), (1-2), (1.7) 
as a function of space x, for various times t ranging from t = to t = 2 (the bottom lines 
correspond to high times). We use the parameters defined in (4.2). 

We then used this basis to compute the reduced solution for a particular (randomly chosen) instance of the 
parameters. The reduced solution was computed in 0.04s, including the time necessary for the online error 
bound computation, shown in Figure 4 (solid line). Our procedure reduces the marginal cost to 15% of the 
original cost, yet providing a certified L 2 relative error of less than l%o. 

Error bound estimation 

Still using the preceding POD basis and instance of the parameters, we compared the online error bound with 
the actual error, for the same parameter set as above. The result is shown in Figure 4. We see that our error 
bound is quite sharp, especially when it follows the decrease in the actual error near t = 1.3. We also checked 



TITLE WILL BE SET BY THE PUBLISHER 



Parameter 


Min. 


Max. 


Parameter 


Min. 
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V 
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1.1 
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.9 


1.2 




1 


1 




.9 
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bi 
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1 


1 


frn 





2 


, J T 
UJ 1 


2 


2 


A i,i 


0.7 


1.3 


W l 


2 


2 







1 


, no 


3 


3 



Table 1. Ranges of the different parameters in the benchmark of Section 4.2. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



0.132 I 1 1 1 1 1 

0.13 ■ 

^ 0.128 ■ \- 

S 0.126 - - 
N 0.124 - 
0.122 ■ 

0.12 1 1 1 1 

0.2 0.4 0.6 0.8 1 

x 

FIGURE 3. POD reduced basis: plots, as functions of space, of the 7 leading POD modes, i.e. 
the Q (i = 1, ... , 7) defined by (2.6), with z t (i = 1, . . . , 7) the leading eigenvectors of M T QM. 
The modes are sorted (from top to bottom, and from left to right) by increasing magnitude of 
eigenvalues. Parameter ranges for snapshot sampling are those in Table 1. 



for the quality of the SCM procedure, by comparing the actual stability constant with the lower bound provl 
by SCM (Figure 5). 
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0.0001 



1e-05 
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1e-07 



1e-08 



1e-09 




L2 error bound (relative) 
L2 actual error (relative) 



0.5 
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9 &k 

Figure 4. Relative L online error bound and actual error. We plot in solid line 

\\u k \ 



and 



in dashed line 



as functions of t — kAt for k = 1, . . . , T. 



9.95 



9.9 



9.85 



9.8 



9.75 



9.7 



9.65 





true stability constant — 




SCM lower bound 






\ 




\ 







0.2 0.4 0.6 0.£ 



1 1.2 1.4 1.6 1.8 2 



Figure 5. True stability constant Cfc denned by (3.3) and SCM lower bound C™ f denned by 
(3.11) as functions of time. We use M = #C = 10 as SCM parameters. 

4.3. Convergence benchmarks 

In order to compare our three basis selection procedures (POD, greedy and POD-Greedy), we have made 
"convergence benchmarks", i.e. representations of the maximal and mean (estimated) error over all timesteps, 
and over a sample of 100 parameters as functions of the size of the reduced basis. The same sample of benchmark 
parameters is used throughout all the procedure. 

Comparison of greedy (with #S = 100) and POD (with S = 60) procedures for TV = 40, T = 2, At = .02, 
n(bo) — n(bi) = ns(f) — nr(f) — n(uo) = 0, with parameters f m and v fixed to unity, and varying uo m S [0, 1] 
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0.1 r 1 1 1 1 

mean relative error on 1 00 values (POD) 

maximal relative error on 1 00 values (POD) 

mean relative error on 1 00 values (GDY) 

\ maximal relative error on 100 values (GDY) 

0.01 - \ 



0.001 - 




1e-06 1 ' ' ' ' 1 

5 6 7 8 9 10 

reduced basis size 

Figure 6. Convergence benchmark 1. We plot (on a logarithmic scale) maximal and mean 
relative online error bounds over a (uniform) random sample of 100 initial values uo m € [0, 1], 
as functions of the reduced basis size N, when the reduced basis is chosen using POD-based 
procedure (POD) or greedy procedure (GDY). Fixed parameters are v = 1 and f m = 1. 



(and thus initial boundary values &om and bi m , moving accordingly to compatibility conditions (1.4)) is shown 
in Figure 6. 

The benchmarking process for greedy took 19.5s of CPU time, the one for POD took 14.98s. The online cost, 
depending only on the size of the reduced basis, is the same regardless of how the reduced basis has been chosen. 
We also see the fast (exponential) convergence of error bound towards zero as N increases. 
Another benchmark was then made, with the same data, except that v = .1 and At = .002. The result is 
visualized in Figure 7. The POD benchmarking process took 349 s of CPU time, the POD-Greedy, 282 s, and 
the Greedy, 470 s. We notice that a smaller viscosity leads to degraded precision of our RB approximation. The 
resulting final basis selected by the POD-Greedy is displayed in Figure 8. 

The POD-Greedy algorithm is run using Pi = 2 and initialized using the full time-discrete trajectory (hence, 
200 vectors) for one parameter value. Using the initial data as initialization for the POD-Greedy algorithm 
and using Pi = 1 may give better overall performance of POD-Greedy (at the expense of an increased offline 
computation time), as the sample of error bounds is updated at every step of the algorithm (instead of every 
other step for Pi = 2). 

4.4. Effect of mesh refinement and penalization constant 

To check for the robustness of our bound, we studied the influence of Af (the number of "full" spatial discretiza- 
tion points) and P on the magnitude and the sharpness of our error bound. We ran the same benchmark as in 
Section 4.2, but with different Af or P. 

We did the test with refined meshes (Af = 200, Af = 800) and we obtained a similar error bound profile; we 
conclude that the sharpness of our error bound is quite insensitive to mesh refinement. 

In Figure 9, we visualize the actual error and the error bound for various values of P (with Af = 60). We see 
that the error bound is tighter and sharper for high values of the penalization constant P. This can easily be 
explained by the fact that, as P — > +oo, the errors at the boundary efc(0) and efc(l) vanish, hence modifying 
all the terms in the error bound (i.e., decreasing rj^, af., fk, £fc\ , G, i5fc, 7fe and increasing Ak) where these 
errors appear. 
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Mean relative error bound (Greedy) 
Mean relative error bound (POD) 
Mean relative error bound (POD-Greedy) 
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Maximal relative error bound (POD) 
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16 



Figure 7. Convergence benchmark 2. We plot (on a logarithmic scale) maximal (bottom) 
and mean (top) online error bounds over a (uniform) random sample of 100 initial values 
uom € [0, 1] , as functions of the reduced basis size N, when reduced basis is chosen using 
POD-based procedure with S = 90, Greedy, or POD-Greedy procedure with Pi = 2. Fixed 
parameters are v = 0.1 and f m = 1. 



4.5. Comparison with existing bound 

We compared our error bound with the bound described in [18], for M = 60, At = .02, T = 2, fixed uq = bo = 
bi = 0, / = 1 and variable v £ [0.1; 1]. The reduced basis is found by POD with S = 90. Figure 10 shows that 
our bound is clearly better than the existing reference bound. 

Besides, we have performed a benchmark over a fixed sample of 100 random values of v uniformly chosen in 
[0.1; 1]. For our bound, the mean error bound is 0.00076, the maximum error bound is 0.02, while for the bound 
of [18], we obtain 0.0041 for the mean bound and 0.25 for the maximum. 

It is easy to see that, when the boundary conditions arc fixed to zero, the recurrence formula for the error bound 
reduces to: 



< 



|e fc _i|| + Atjjrfcj 
1 + C7 fe At 
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FIGURE 8. Reduced basis selected at the final step of the POD-Greedy procedure carried out 
for the previous benchmark (Figure 7). Basis elements are plotted as functions of space, and 
are sorted (from top to bottom, left to right) in order of selection in the greedy procedure. 



while the error bound described in [18] has the following expression: 
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L2 error bound (relative) 
L2 actual error (relative) 



L2 error bound (relative) 
L2 actual error (relative) 



L2 error bound (relative) 
L2 actual error (relative) 




L2 error bound (relative) 
L2 actual error (relative) 



2 £ k 

Figure 9. Relative L online error bound and actual error. We plot in solid line , , and 



in dashed line 



\u k — u k \ 



as functions of t = kAt for k = 1, . . . , T, and Af — 60 and, from left 



to right and top to bottom, P = 10 2 , P = 10 4 , P = 10 5 and P = 10 



12 



where the modified stability constant Ck reads: 



C k = inf 

v£X 



4c(u h , v, v) + ua(v, v) 



We recall that our stability constant Ck is given by: 



Ck = inf [2c(u fc , v, v) + i/a(v, v)] 
uex ,|M|=i 



The ^-dependence of the bound can explain why our bound is better, especially for small values of v. Besides, 
the derivation of our error bound makes lesser use of inequalities (for instance, we do not make use of Young's 
inequality at the beginning of the proof, each inequality used is a potential source of optimality loss) and keeps 
treating more terms. 
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0.0001 



1e-05 

0.5 1 1.5 2 

t 

Figure 10. Comparison with the existing error bound for reduced basis Burgers equation 
(reference bound). We plot the actual error, the existing (reference) error bound and our error 
bound as functions of time. We took v = 0.1. 

Conclusion 

We have presented a certified procedure for low marginal cost approximate resolution of the viscous Burgers 
equation with parametrized viscosity, as well as initial and boundary value data. This procedure makes use of 
a reduced basis offline/online procedure for a penalized weak formulation, an efficiently computed error bound 
in natural L 2 norm (made possible by the successive constraints method (SCM)), and three procedures at hand 
for choosing a basis to expand reduced solutions in. 

Our procedure becomes less useful when the ratio time/ viscosity increases, as this degradates the stability 
constant Ck- Another limitation of our method, for one willing to use it for large times, is that the online 
procedure complexity still depends on the temporal discretization step. However, our numerical experiments 
show a substantial decrease in marginal cost when using reduced basis approximation, as well as efficiency (both 
in terms of sharpness and computation time) of the provided error bound for moderate viscosities. This decrease 
in the cost is made possible by the fact that online procedure has a complexity that is independent from the 
number of spatial discretization points. 

Appendix A. Proof of Theorem 3.1 

Proof. Subtracting left-hand side of (2.2b) from both sides of relation (3.1) yields that for k = 1, . . . , T, the 
error at time tk- &k — — u satisfies, for every v € Xq: 

-^((e k ,v) - (e k -i,v)) + c{u k e ,ul,v) - c{u k ,u k ,v) + va(e k ,v) = r k (v) (A.l) 

We write: 

e k =e fc (O)0o + e fc (l)(/w + ei (A.2) 
with e% e X Q = {v e X st. v(0) = v(l) = 0}. And, by applying (A.l) with v = e%: 

4r((efe, el) - (e k -i, e z k )) + c(u k e ,u k e , e z k ) - c{u k , u k , e z k ) + va{e k ,e z k ) = r k {e z k ) (A.3) 
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Since: 

(e k ,e z k ) = e fe (O)(0 o ,eI) + e k (l)(^,e z k ) + \\e z k \\ 2 

= ||e£|| 2 + e k (0)el (±A (<f> , fa) + e k (l)e% (l - i ) 



and that, for every v <G X: 



c(ul,u k e ,v) - c{u\u\v) = -\ f ((„*)» - (^) 2 ) £ 



9x 



1 Z" 1 _ fe \ <9t> 

= -27 ^ + ^ + ^ 

= 2c(u k ,e k ,v) + c(e k ,e k ,v) 

We have that (A. 3) implies: 

i-(||e^|| 2 -(e fc _ 1 ,e^))+^(e fc , e z fc ) + ^ [e k {Q)el (j^j (fa, fa) + e k (l)e% (l-jf) (fo,4>tf-i) 

= r k (e z k )-c(e k ,e k ,e%) (A.4) 

We are now willing to find a lower bound for the left-hand side of (A.4) and an upper bound for its right-hand 
side. 

Lower bound for LHS. From Cauchy-Schwarz inequality: 

-(e k -i,et) > -||e fc _i||||e||| 

and, by the triangle inequality: 

||e fe ||-r ?fe <||eI||<||e fe ||+ % (A.5) 

because -q k = |e fc (0)| ||</> || + |e fc (l)| \\fav\\. 
So: 

-(e k -i,e z k ) > -||e fc _i||(||e fc ||+» ft ) (A.6) 

We also have: 

II4II 2 > Ik-f + \\e k e\f - 2 ||e fe || \\e k e* k \\ > \\e k f 2 ||e*|| m (A.7) 
because of Cauchy-Schwarz and triangle inequalities. 
Besides, 

0<||ejE|| < ||e^-e fe || + ||e fc || < % + ||e fe || 

so that: 

||4|| 2 <(|| efc ||+%) 2 (A.8) 

Using the bilinearity of ip k , we have: 

ip k {e k , el) = ip k (e z k , e z k ) + e fc (O)^ fc (0 o , <h)ef. (J^J + e k (l)ip k (fas , <t>N-i)e k ^1 - (A.9) 

because ip k (4>o,<f>j\f) = ^ki&Afi&o) — 0, since <f>Q and 4>M have no common support, ip k (4>Q,4>j) = for j > 1, 
and i/j k (fav, fa) = for j < Af - 1. 
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From the definition of the stability constant C k : 



So that, thanks to (A. 7) and (A. 8), 

ipk(e z k ,e z k ) > 

That is 



M4,eZ)>C k \\e' k \\' 



Cfc||e fc |r- 277fcCfc||e fc || 



if C k > 



C k \\e k \\ 2 + 2 Vk C k \\e k \\ + C k rg if C k < 
Me%,e%) > C k \\e k \\ 2 - a k \\e k \\ - [C k }_ 4 



from the definition of o~ k . 
We have: 



<£\\e z k \\<£\\e k \\+£r) k 



and by symmetry: 

el (l-^=) <£\\e z k \\<£\\e k \\+£r ]k 
so that, combining (A. 11), (A. 12) and introducing f k : 



e fc (O)V'fc(0o,^i)eI ( ) +e fc (l)Vfc(<?W,<?W-i)eI ( 1 - 



< ||ejfe|| fk + Vkfk 



Now we can say that: 



AT 



e/c(O)-0fc(0o,^i)e| 77 + e k (l)ip k (<f)j^, 4>^-i)e k 1 



Af 



> - ||efe|| /fe - r?fc/fe 



And we also have: 



e fc(°) e fe ( JJ ) (0o,0i) +e fc (l)el, ( 1- — f ) {<j>M-i,4>M) 



(A.10) 

(A.11) 
(A.12) 



(A.13) 



< |e fc (0)|(£||e fc || +£77 fc )(0o,<£i) 

+ |e fc (l)|(f ||e fe || +£ Vk )(<t>M,<t>M-i) 
= £(\e k (0)\ (to, + |e fc (l)| («W, <^-i» ||e fc || 

+ ^ fe (|e fe (0)| (0o, ^i) + |e fe (l)| (0at,^V-i)) 
= f(0o,0i)(| efc (O)| + |e fc (l)|)||e fc || 

+ ^ fc (to,^i)(|e fc (0)| + |e fe (l)|) (A.14) 

since (<fo,0i) = {4>mAm-i) by symmetry. 

Thus, thanks to (A.7), (A.6), (A.9), (A.10), (A.13) and (A.14), the left-hand side of (A.4) is greater than: 



s +<fc)iie„ir- 



2 /2 % + ||e fc _i|| +£(0 o ,0i)(|e fc (O)| + |e*(l)|) 



Af 



+ <7 fc + / fc ||e fc || 



||e fc -i|| +g^fc<0o,0i)(|e fc (O)| + |e fc (l)|) 
At 



[C k ]_vl-Vkfk (A.15) 



Upper bound for RHS. We have: 

c(e k , e k , e%) = c(e k ,e k ,e k ) - e fc (0)c(e fc , e k , <j> ) - e k (l)c(e k ,e k , 4>m) 
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but: 



So: 



c(efe, ek, e k ) 



1 f 1 2 de k 
°' k dx 
1 /^[(e,) 3 ] 



6 ./n 



-i(( efc (l)) 3 ~(e fe (0)) 3 



c(e fe! e h , ej) = - i ((e fc (l)) 3 - (e,(0)) 3 ) - (e fc (0) £ e 2 fc ^ + e fe (l) jf* e 2 ^) 
= -6((e fc (l)) 3 -(e fe (0)) 3 )+^(e fe (0) ^ e 2 , - e fc (l) ^ 



L-l/JV 



Since, for all x £ 



we have, thanks to 



6k[ Af 



e k (x) = e fc (0) +7V ( e k 
<£||e*||: 



ejfe(O) a; 



1/JV 

Afxe k (0)l e 2 



e k (0)\ 



|e fe (0)| 



e fc (^) 2 + e fc (0) 2 , 2|e fc (^)e fcl 



Pll III ""i , g2 II^H 2 , e fc(°) 2 , 2 Ie fc (0)| ^ 
< |e fc (0)| [£ \\e k \\ \e k (0)\ H 1 1 £ \\e k \ 



3 3 3 

3 



^ g 2 |e fc (0)| I, ,,2,5 a M , |e fc (0)| 
< 3 l|efc|| +-|e fc (0)| t\\e k \\ + — 



As a similar computation can be worked out for 



Mxe k {l) 



l-l/JV 



el 



we have: 



e^O) I el^+e k (l)l ) < \\e k \[' + \\e k \\ + £ 



dx 



1 2 d ^\ 

et— — 



cA ||„ II 2 . C B 



dx I 



We also have, thanks to (A. 5): 



where: 



|r*(e£)| < ||r fe || ||e||| < ||r fc || ||e fc || + ||r fc || m 



Ikfello = su p r k(v) 

uex ,||«||=i 



Hence, the right-hand side of (A. 4) is less than: 



^ \\ek\\ 2 + (\\r k \\ + \\e k \\ + - \e k (lf e k (0f\ + % + \\r k \\ r, k 
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Conclusion. Now (A. 4) implies, thanks to (A. 15), and (A. 16): 

A||e fc || 2 -£ fc ||e fe || - 7fc <0 (A.17) 

Viewing left-hand side of (A.17) as a (convex, thanks to our hypothesis (3.4)) quadratic function Q of ||efc||, 
whose discriminant is T>k, equation (A.17) implies that, if T>k > 0, ||efe|| is smaller than the greatest real root of 
Q, that is: 

B k + ^JW k 

If T>k < 0, then necessarily j k < (as Ak is positive). Hence (A.17) implies: 

A||efc|| 2 -£fc||efc|| <0 



and so: 



A 



k 
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